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ABSTRACT 

Coronal holes are the source regions of the fast solar wind, which fills most 
of the solar system volume near the cycle minimum. Removing stray light from 
extreme ultraviolet (EUV) images of the Sun's corona is of high astrophysical 
importance, as it is required to make meaningful determinations of temperatures 
and densities of coronal holes. EUV images tend to be dominated by the com- 
ponent of the stray light due to the long-range scatter caused by microroughness 
of telescope mirror surfaces, and this component has proven very difficult to 
measure in pre-flight characterization. In-flight characterization heretofore has 
proven elusive due to the fact that the detected image is simultaneously nonlin- 
ear in two unknown functions: the stray light pattern and the true image which 
would be seen by an ideal telescope. Using a constrained blind deconvolution 
technique that takes advantage of known zeros in the true image provided by a 
fortuitous lunar transit, we have removed the stray light from solar images seen 
by the EUVI instrument on STEREO-B in all four filter bands (171, 195, 284, 
and 304 A). Uncertainty measures of the stray light corrected images, which in- 
clude the systematic error due to misestimation of the scatter, are provided. It is 
shown that in EUVI, stray light contributes up to 70% of the emission in coronal 
holes seen on the solar disk, which has dramatic consequences for diagnostics of 
temperature and density and therefore estimates of key plasma parameters such 
as the plasma (3 and ion-electron collision rates. 
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INTRODUCTION 



One of the longest standing puzzles in astrophysics concerns the processes that energize 
the solar wind. The fast solar wind, which has a speed of ~ 800 km/s near the Earth, 
emanates frorn EUV and X-ray faint regions in the Sun's atmosphere called coronal holes 
( iKrieger et al.lll973l ). Physics-based modeling efforts to identify the processes t hat heat and 
accelerate the solar wind have met with very limited success (jCranmerl |2010| ) . It is likely 
that much more observational input will be required to adequately constrain the numerical 
models used to investigate this question. NASA's EUV imaging instruments (SOHO/EIT, 
TRACE, STEREO/EUVI, and SDO/AIA) provide a comprehensive data set covering 11/2 
solar c ycles that can prov ide powerful constraints on temperatures and densities in coronal 
holes (IFrazin et al.ll2009l ). However, coronal holes and other faint structures are severely 
contaminated with stray light, which must be corrected before these vast data sets may be 
utilized for this purpose. 

All EUV imaging instruments are afflicted to some degree by long-range scattering, 
which distributes a haze of stray light over the whole imaging plane. This haze is not no- 
ticeable in the brighter areas of the image, but can completely overwhelm the emissions 
of faint regions. Stray light arises from two sources: non-specular reflection from micro- 
rough mirror surfac es, and diffraction due to the pupfl and a mesh obstructing the pupil 
( iHoward et al.ll2008l ). Here we describe a method for, and results of, correction of stray light 
in STEREO-B/EUVI, henceforth abbreviated as EUVI-B. 

Stray light correction requires determination of the instrument's point spread function 
(PSF). The four EUVI channels (171, 195, 284, and 304 A) have different optical paths and 
therefore different PSFs. Except for pixel-scale variations due to optical aberration (whereas 
the scattering of interes t here is significant o ver ranges of hundreds of pixels), the EUVI PSFs 
are spatially invariant (IHoward et al.ll2008l ). Thus, stray light contami nation is modeled b y 
convolution with the PSF, and the correction process is deconvolution (jStarck et al.ll2002[ ). 



The instrument PSFs are difflcult to characterize experimentally due to the lack of a 
sufflciently strong EUV source, so they must be determined primarily from in-flight obser- 
vations. So lar flare images caii provide signiflcant information about the entrance aperture 
diffraction ( iGburek et al.ll2006l ). as diffraction orders are easily visible around the flare. Un- 
fortunately, most of the scatter is too diffuse to be observed clearly around flares. The best 
information about diffuse scatter is obtained from transiting bodies that do not emit in the 



EUV, so any apparent emission is instrumental in origin. DeForest et al. (IDeForest et al. 



20091 ) used a Venus transit to estimate the mirror scattering in the TRACE instrument by 
fitting a truncated Lorentzian. Here, we present the first self-consistent determination of the 
mirror scattering via blind deconvolution, in which both the true solar image and various 
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PSF parameters are taken to be unknown simultaneously. The information making this effort 
successful comes from calibration rolls and the Feb. 25, 2007 STEREO-B lunar transit, each 
exposure of which provided about 50,000 lunar disk pixels illuminated only by instrumental 
effects. 



2. BLIND DECONVOLUTION METHOD 

We characterize scattering in a given EUVI filter band with a three-component PSF. The 
first two components, /i^ and /i^, account for diffraction through the pupil and Pupil mesh, 
and were determined analytically using Fraunhofer diffraction theory f Goodmanlll996l ) . The 
third component, /i™, accounts for the remaining diffuse scatter, much of which derives from 
the mirror microroughness. The formula for this component contains free parameters 
which we determine empirically by blind deconvolution, and we write h"^ to indicate the 
dependence. The total PSF is the convolution of these three components: 

h^ = ¥*hP* /i™, (1) 

where 

{u* v){x) = ''^u{x — x')v{x') (2) 

x'£l 

denotes the convolution of two discrete functions u and v over the index set / of the 
2048 X 2048 array of CCD pixels. (We set u{x - x') = at pixels x - x' ^ I.) This 
model assumes that each component is independent, neglecting phase correlations that arise 
as light propagates from the entrance aperture to the primary mirror. To determine the 
proportionality constants of the three PSFs, we assume that each sums to unity. 

We used the EUV mirror literature to help us choose an appropriate parametric model 
for the empirical PSF component h'^. In a typical EUV scatter model, a si gnificant fraction 



of th e light is not scattered, while the rest is scattered by broad wings (IKrautschik et al. 



2OO2I ). Up to scaling and normalization constants, these wings are described by the power 



spectral density (PSD) of the mirror sur face height function, and this PSD has been directly 



measured for mirrors similar to EUVI's (IMartmez-Galarce et al.ll2010l ). A log-log plot of the 
measured PSD versus spatial frequency is roughly piecewise linear, implying that is a 
piecewise power law whose exponent depends on pixel distance r from the origin. 

Accordingly, a parametric formula for a family of piecewise power laws was used to 
describe h^. A series of breakpoints 1 = tq < ri < r2 < . . . < r?, < 00 was chosen, and the 
number of breakpoints, b = 8, was selected according to a goodness of fit criterion (we 
performed several fits using different values of b, and the fit did not improve significantly 
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for b > 8). On each subinterval [rj_i,rj), the formula is given by Pa,^{r) oc r~^\ where 
Pi > 0, and /3 = . . . , Pb). The parameter a represents the fraction of non-scattered hght: 
Pa,/3(0) = a. This radially symmetric profile did not adequately describe the anisotropic 
scatter observed in the transit and calibration roll images, even after removal of the scatter 
due to and h^. We therefore generalized the formula to allow to have an elliptical 
cross section. Let Mg^e denote the 2x2 matrix that dilates the plane by a factor of s along 
a line rotated 6 radians counterclockwise from the horizontal axis; then 



h':^{x)=pa,(i{\\Ms,e-x\\2 + 6), (3) 

where the constant S = 1 was added to avoid the power law's singularity at the origin. The 
free parameters of the PSF model are then if = {(3, a, 6, s) G M*', where p = b + 3 = 11. 

To determine the PSF from the lunar transit data, we described the observation process 
by a statistical image formation model and sought a maximum likelihood estimate of the 
PSF under this model. We let u^^^'^ denote the ideal image; that is, u^^^'^{x) is the expected 
solar photon count that would be measured at pixel x by an ideal instrument. Due to scatter 
and the Poisson photon arrival process, the actual number of photon arrivals is a Poisson 
random variable with expected value h^^^° * The difference between the expected and 

the actual photon count is the photon noise riphot- The observed photon count / deviates 
from /i*''™ * ytruc _|_ j2pjjQ^ due to CCD dark current and read noise need- Based on histograms 
of dark images, we find need is reasonably modeled as a Gaussian white noise process with 
standard deviation (Teed ~ 1 digital number (DN). The variance must be divided by the 
photometric gain factor (the recorded DN per incident photon) to obtain the CCD noise 
level in units of photons. Combining nphot and rieed into a single variable n, we obtain the 
statistical image formation model 

/ = h'''''' * u*''"" + n. (4) 

We would like to find a maximum likelihood estimate of the PSF under the model of Eq. 
but the Poisson- Gaussian distribution of n leads to a difficult large-scale nonlinear opti- 
mization problem. We therefore apply a variance stabilizing transform (VST) to make n 
approximately standard normal (mean zero and variance one), which leads to an easier 
least-squares problem. A VS T for a Poisson- Gauss ian random variable where the Gaussian 



has variance is derived in ( iMurtagh et al.lll995l ). Their full formula allows for CCD bias 



and non-unity camera gain; however, we correct for these in pre-processing, so the formula 



cr^ = we obtain 



simplifies to vst(/) = 2^/J + | + ^2. if the VST is applied to both sides of Eq. with 

vst(/) ^ vst(/i*™'= * Mt™<=) + nvst, (5) 
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where, for each x, n^stix) is roughly standard normal (iMurtagh et al.lll995l ). The approxi- 
mation begins to break down when / < 5 photons, but this only occurs far from the solar 
disk in the lunar transit images. In these low-intensity regions we replace the observed image 
pixel value with a local average of sufficient size that the total photon counts exceed 5 in the 
averaged neighborhood. 

We found an estimate of the EUVI-B PSF h^^^"^ by solving a blind deconvolution 
problem on a series of 8 images /i,...,/8 from the lunar transit series. We assume the 
scatter- free images m*"""® are positive everywhere and zero on the lunar disk pixels Zi, leading 
to the constraints 

nf->0, nf-(Z,) = 0, z = l,...,8. (6) 

The lunar disk was identified by detecting its edge pixels with gradient thresholding, then 
fitting a circle to the detected edge pixels. An approximate maximum likelihood estimate 
for the PSF and images {ui} is obtained by solving 

8 

minimize \ \\vst{h^ * Ui) — vst{fi)\\'^ 

(7) 

subject to Ui > 0, Ui{Zi) = for all i, 
f3>0, < a < 1. 

Since each Ui is a 2048 x 2048 image, this problem has over 32 million variables and is in- 
tractable by general-purpose numerical meth ods. We solved this pro blem with a customized 



variant of the variable projection method fiGolub Sz Pereyral |2002| ). which is efficient for 
problems of this kind. 

Fig. [1] shows a log- log plot of the 171 A profile function obtained by blind deconvolu- 
tion, along with views of the total PSF and its core. Note that the profile is constant for 
200 < r < 1200 pixels, implying considerably more long range scatter than a constant power 
law exponent would predict. The total PSF is dominated by the elliptical power law decay, 
although the diffraction orders from the pupil and mesh PSFs are visible near the PSF core. 

Once we determined the PSF from blind deconvolution, we were able to correct any 
observed EUVI-B image / via more standard deconvolution methods. Let C{h) : M.^ M.^ 
denote the linear operator that convolves an input image M with a fixed PSF /i: C{h)u = h*u. 
The stray light correction of a given image / is obtained by applying the inverse of C{h) to 
/: 

u = C{h)-'f. (8) 

The inverse operator does not exist for arbitrary h. However, our PSFs have an origin value 
h{0) > 1/2, which makes C{h) diagonally dominant and invertible. The inversion of Eq. ([8]) 
can be performed with the conjugate gradient method in less than a minute on a laptop. 
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3. MODEL VALIDATION AND ERROR ANALYSIS 

A first test of our stray light correction was performed by applying cross-validation (CV) 
to the lunar transit series. For each of the 8 images fi, we determined a PSF hi by solving 
Eq. ([7]) with fi removed from the dataset. We then calculated the stray light corrected CV 
image u* = C{hi)~^ fi. Note that both hi and u* are calculated without any assumption of 
a zero lunar disk in image so the values of u* on the lunar disk represent an independent 
check on the correction's effectiveness. In Fig. [2] we compare the lunar disks of fi and u* 
in 171 A. We find that the lunar disk values after deconvolution are very strongly clustered 
near zero, as seen in a histogram of the lunar disk values {bottom right). 

Calibration roll images provide direct evidence of anisotropic scatter and our deconvo- 
lution's ability to correct it. We will present the analysis for 171 A here, leaving the other 
bands to future publications. On Nov. 8 2011, STEREO-B executed a 360° calibration roll, 
and in each band, EUVI-B acquired 9 solar images at roll angles of 0, 60, 90, 120, 180, 240, 
270, 300, and 360° relative to the pre-roll position. These images are useful because the 
direction of anisotropic scatter rotates with STEREO-B, introducing discrepancies between 
the images. To pick out these discrepancies, an 8 x 8 boxcar was applied to reduce noise 
and the images were rotated into a common solar coordinate system. This procedure was 
repeated on stray light-corrected roll images, resulting in series {/p}p=o ^^'^ {'"p}p=o indexed 
by roll angle p. A haze of scattered light can be observed off the limb in the the image /o 
(Fig. m top left), which is much reduced in the corrected image uq [middle left). 

Anisotropic scatter rotation was tracked using the difference images Afp = fp — fo and 
Aup = Up — Uq. To avoid analyzing regions of the Sun that changed significantly over the 
course of the roll, we examined the difference A/360 between the pre- and post-roll images, 
and all pixels where IA/360I > 1 DN were masked out of the difference images. In A/90 
[top right), we observe a 'dark axis' and a ~ 90°-rotated 'light axis' corresponding to the 
preferential scattering directions for /q and /90 respectively. (These regions are very diffuse, 
so the separation angle can only be estimated up to about ±10°.) These axes are eliminated 
in Augo because the stray light correction greatly reduces the anisotropic haze. Similar 
reductions are observed at other p values. 

To estimate the error in the corrected images, we decompose the total error e = u — u^^^*^ 
into components due to noise and PSF error, and estimate their contributions separately. 
To obtain the decomposition we substitute Eq. (jlj) into Eq. and introduce the PSF error 
variable 6h = h — /i*''"®, obtaining 

e = u- M*™" = C(/i)"^C(5/i)m*™'= + C{h)~^n. (9) 

The right hand side's first term gives the error due to PSF misestimation, which we name 
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^psf. The second, enoise, is the noise in the corrected image. 

To estimate ^psf, we first apply an 8 x 8 moving average filter to u and / to reduce noise, 
so that u — u^^^^ ^ Epsf. (Averaging also reduces ^psf on small scales, potentially leading 
to underestimation of the error. However, we are only interested in the large-scale effect of 
stray light, and on large scales, epsf is unaffected by averaging.) Next, we assume that 5psf 
can be bounded by a quantity (Tpsf which is proportional to the magnitude of the stray light 
correction \u — f\: 

\u-u''^''\<B\u-f\ = ap,i. (10) 

We use the bounding term B\u — f\ because it shares two distinctive properties with 5psf: 
both are approximately invariant if u^^^'^ is changed by scaling or addition of a constant. 
These invariance properties help ensure that an empirical bound based on lunar transit 
images alone will continue to be valid for other solar images. 

We determined the constant B by requiring that Eq. ffTOl) hold true over almost all of 
the lunar disk pixels of fi and u*. Setting u = u*, f = fi, dividing both sides by \u* — fi\, 
and noting that uf^"^ = on the lunar disk, we find \u*\/\u* — fi\ < B is required on the 
lunar disk pixels. We call the left-hand side ratio b*, and collect the values of b* on the 8 
lunar disks into a histogram. A normalized histogram for 171 A is given in Fig. [3l right, 
and values at the 68*'^, 95*^, and 99.7*^ percentiles, corresponding to the first three standard 
deviations of a Gaussian, are 0.08, 0.13, and 0.16. We set B equal to the 95*'* percentile of 
the histogram, which gives a conservative 2a error bound robust to outliers. To estimate the 
contribution of noise to the error, we analytically calculate its covariance matrix using the 
known variance of n. The diagonal of this covariance matrix, cr^oisc^ is our estimate of the 
squared error due to noise. We add the two independent errors in quadrature to obtain the 
final error estimate a: = cr^„f + crf-^. 



4. RESULTS 

Here we note just a few of the effects of stray light correction on faint regions. Fig. O 
shows a coronal hole before and after correction. The stray light corrected coronal hole is 
significantly dimmer: the percent change {u — f)/f ranges from —40 to —70% over most of 
the coronal hole, implying that 40-70% of the apparent coronal hole emissions are stray light. 
The plot of intensity versus position (right) contains the error bars a. Note that they are 
small relative to the size of the correction, which is representative of typical results. Stray 
light correction also removes the haze seen above the limb (Fig. HI top and center left), often 
reducing the intensity by 75% or more (bottom left). A filament just above the solar center 
is also reduced in intensity by 50-70% after stray light correction (bottom left). 
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These large downward corrections to faint regions have major impacts on the plasma 
diagnostics available from EUV images, which in turn are related to the electron temperature 
Te and density rtf.. Assuming electron impact excitation, the intensity of the emission in 
images from the k^^ spectral band is Ik oc Jj^^g d/ ng(/)(5fc(Te(/)), where Qk{Tg) comes from 



an op tically thin p 



sight (IBrown et al 



asma emission model and /lqs represents integration along the line-of- 



19911). This relationship may be utilized to create 3D tomographic maps 



of He and Tg ( iFrazin et al.l 120091 ) . A downward correction of the observed intensity causes a 
proportional reduction of the estimated value of nl, indeed, analysis of the 171, 195 and 284 
A intensities show a downward revision of the coronal hole column density [Ji^Qa^dln'Kl)]^^'^ 
of ~ 40%. The effect on the estimated Tg is more complex. The removal of stray light 
from the off-limb causes a dramatic steepening of the profile function ne{h) (where h is 
the height above the photosphere). The specific impact of stray light correction, including 
constraints on solar wind models from the corrected profiles ne{h), Te{h), is currently under 
investigation. Obvious consequences include reduction of the plasma /3, electron-ion collision 
rates, and the mass of solar wind plasma requiring acceleration. 



5. CONCLUSION 

We have obtained PSFs for all 4 bands of STEREO-B/EUVI using lunar transit data, 
which enables us to correct all EUVI-B images for stray light and provide uncertainties. 
Similar methods may be applied to treat the stray light problems in the other solar EUV 
imaging instruments (SOHO/EIT, TRACE, STEREO-A/EUVI, SDO/AIA). This work and 
its heliophysical implications will be reported in more detail in upcoming publications. 

We thank Jean-Pierre Wuelser for continuous assistance with EUVI technical issues; 
Frederic Auchere and Raymond Mercier, for useful discussions of PSF modeling; and Simon 
Plunkett, for granting our request for EUVI calibration roll data. 
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Fig. 1. — An overview of the 171 A PSF structure. All distances are in pixels. Left: log-log 
plot of power- law profile function Pa,i3{f) with breakpoints marked. Center: total PSF 
with elliptical contour highlighted in white (logarithmic color scale). Right: PSF core. The 
centrally emanating streaks are diffraction orders from and . 
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Fig. 2. — The lunar transit before and after stray light correction (171 A). The lunar and 
solar limbs are outlined for reference. Top left: The lunar disks from the eight transit 
images /j superimposed on the first image of the series, /i. The colorbar is in photons 
per second (ph/sec) and has a low upper limit to show the stray light on the lunar disks. 
The black line segment identifies a series of pixels whose intensity values are plotted to the 
right. Bottom left: The lunar disks from the CV images u* superimposed on Ui. Top right: 
Intensity along the vertical black line segment before correction (squares) and after (solid). 
Bottom right: Normalized histograms giving the distribution of lunar disk intensities before 
correction (squares) and after (solid). 
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Fig. 3. — Left: Map of the b* values on the 8 lunar disks. Right: Normalized histogram 
giving the distribution of the b* values. The dashed vertical lines indicate the 68*^ , 95*'', and 
99. T*'^ percentile values. 
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Fig. 4. — Analysis of calibration roll images before and after stray light correction. Left: The 
first roll image before stray light correction (top), after correction (middle), and fractional 
change ((mq — /o)//o) bottom). The colorbars for /o and Uq are in units of log;^o DN. Right: 
Sample difference image before correction (A/90, top) and after (Awgo, middle), units of DN. 
The two black lines denote the axes of preferential scatter for /o and /go- 
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Fig. 5. — Stray light correction of an on-disk coronal hole observed by EUVI-B on Nov. 21, 
2008. Left and center: The observed data / and corrected data u (ph/sec). Right: Line 
plots of / (crosses) and u (solid line) along the white line segment. The horizontal axis is in 
pixels relative to the segment's left endpoint. Error bars are given every 20 pixels. 



